standard_theme <- theme_bw() +
theme(plot.title = element_text(size = 30, hjust = 0.5),
plot.subtitle = element_text(size = 25, hjust = 0.5),
plot.caption = element_text(size = 20),
axis.title = element_text(size = 25),
axis.text.y = element_text(size = 20),
axis.text.x = element_text(size = 20),
legend.title = element_text(size = 20),
legend.text = element_text(size = 22),
strip.text = element_text(size = 22))
theme_set(standard_theme)
set_flextable_defaults(table.layout = "autofit")
stat <- list(
"Mean" = ~mean(., na.rm = TRUE),
"Std_dev" = ~sd(., na.rm = TRUE)
)
stat_ph <- list(
"Среднее" = ~mean(., na.rm = TRUE) %>% round(4),
"Медиана" = ~median(., na.rm = TRUE)%>% round(4),
"Стандартное отклонение" = ~sd(., na.rm = TRUE) %>% round(4),
"Коэффициент вариации (CV), %" = ~round((sd(., na.rm = TRUE))/(mean(., na.rm = TRUE)) * 100, 2)
)
ethanol <- read_delim("./data/raw/Jones-1984.csv")
## Rows: 480 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (5): study_id, route, food, site, sex
## dbl (4): subject_id, time_min, conc_g_per_L, dose_g_per_kg
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
ethanol <- ethanol %>%
mutate(across(where(is.character), as.factor))
summary(ethanol)
## study_id subject_id time_min conc_g_per_L dose_g_per_kg
## Jones_1984:480 Min. : 1.00 Min. : 0 Min. :0.0000 Min. :0.68
## 1st Qu.:12.75 1st Qu.: 60 1st Qu.:0.2575 1st Qu.:0.68
## Median :24.50 Median :150 Median :0.5300 Median :0.68
## Mean :24.50 Mean :180 Mean :0.4966 Mean :0.68
## 3rd Qu.:36.25 3rd Qu.:300 3rd Qu.:0.7300 3rd Qu.:0.68
## Max. :48.00 Max. :420 Max. :1.3000 Max. :0.68
## route food site sex
## PO:480 fasted:480 capillary:480 male:480
##
##
##
##
##
str(ethanol)
## tibble [480 × 9] (S3: tbl_df/tbl/data.frame)
## $ study_id : Factor w/ 1 level "Jones_1984": 1 1 1 1 1 1 1 1 1 1 ...
## $ subject_id : num [1:480] 1 1 1 1 1 1 1 1 1 1 ...
## $ time_min : num [1:480] 0 30 60 90 120 180 240 300 360 420 ...
## $ conc_g_per_L : num [1:480] 0 0.94 0.81 0.73 0.65 0.56 0.44 0.34 0.17 0.01 ...
## $ dose_g_per_kg: num [1:480] 0.68 0.68 0.68 0.68 0.68 0.68 0.68 0.68 0.68 0.68 ...
## $ route : Factor w/ 1 level "PO": 1 1 1 1 1 1 1 1 1 1 ...
## $ food : Factor w/ 1 level "fasted": 1 1 1 1 1 1 1 1 1 1 ...
## $ site : Factor w/ 1 level "capillary": 1 1 1 1 1 1 1 1 1 1 ...
## $ sex : Factor w/ 1 level "male": 1 1 1 1 1 1 1 1 1 1 ...
ggplot(ethanol, aes(x = time_min, y = conc_g_per_L)) +
geom_line(aes(group = subject_id)) +
geom_smooth(method = "gam") +
scale_y_continuous(breaks = seq(0, 1.4, 0.1)) +
scale_x_continuous(breaks = seq(0, 420, 30)) +
# coord_cartesian(xlim = c(60, 420)) +
labs(title = "Спагетти-плот - кривые концентрации время",
subtitle = "Для демонстрации тренда применён метод GAM",
caption = "Jones-1984",
x = "Время, мин",
y = "Концентрация, г/л")
## `geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'
## Общая кривая (среднее плюс стандартное отклонение)
ethanol %>%
mutate(time_min = time_min %>% as.character) %>%
group_by(time_min) %>%
summarise(across(conc_g_per_L, stat, .names = "{.fn}")) %>%
ungroup() %>%
mutate(time_min = time_min %>% as.numeric) %>%
ggplot(aes(y = Mean, x = time_min)) +
geom_point() +
geom_line(aes(group = "")) +
geom_errorbar(aes(ymin = Mean - Std_dev,
ymax = Mean + Std_dev),
width = 20) +
scale_y_continuous(breaks = seq(0, 1.4, 0.1)) +
scale_x_continuous(breaks = seq(0, 420, 30)) +
labs(title = "Среднее ± стандартное отклонение\nдля каждой временной точки",
caption = "Jones-1984",
x = "Время, мин",
y = "Концентрация, г/л")
ggplot(ethanol, aes(x = time_min, y = conc_g_per_L)) +
geom_smooth(method = "lm", aes(group = subject_id), se = FALSE) +
scale_y_continuous(breaks = seq(0, 1.4, 0.1)) +
scale_x_continuous(breaks = seq(0, 420, 30)) +
labs(title = "Линейная регрессия концентрации для каждого субъекта",,
caption = "Jones-1984",
x = "Время, мин",
y = "Концентрация, г/л")
## `geom_smooth()` using formula = 'y ~ x'
#k - это насколько мы прибавляем к индексу tmax
pharma_fun <- function(conc, t, dose, k = 2) {
# Моделирование для подсчёта beta, C0, Vd
t_sample <- {{t}}[(first(which.max({{conc}})) + k):length({{t}})]
conc_sample <- conc[(first(which.max({{conc}})) + k):length({{conc}})]
dose <- first({{dose}})
model <- lm({{conc_sample}} ~ {{t_sample}})
Beta <- -as.numeric(model$coefficients[2])
C_0 <- as.numeric(model$coefficients[1])
V_d <- dose / C_0
# Подсчёт AUC0-t
# if (length({{t}}) < 2) AUC <- NA
# if (anyNA({{t}}) | anyNA({{conc}})) AUC <- NaN
t0 = {{t}}[-length({{t}})]
t1 = {{t}}[-1]
c0 = {{conc}}[-length({{conc}})]
c1 = {{conc}}[-1]
AUC = sum(((c0 + c1)/2)*(t1 - t0))
# Подсчёт Cmax
cmax <- max(conc)
# Подсчёт Tmax
tmax <- t[first(which.max(conc))]
return(list(Beta = Beta,
C_0 = C_0,
V_d = V_d,
AUC = AUC,
Cmax = cmax,
Tmax = tmax))
}
ethanol %>%
arrange(time_min) %>%
nest(.by = subject_id, .key = "subject_data") %>%
mutate(pharma = purrr::map(subject_data,
~pharma_fun(.$conc_g_per_L,
.$time_min,
.$dose_g_per_kg))) %>%
unnest_wider(col = pharma) %>%
select(-subject_data) -> pharmacokinetics_data
pharmacokinetics_data <- pharmacokinetics_data %>%
mutate(Beta = Beta * 60) %>% # перевожу в г/л в час
rename(`β, г/л в час` = Beta,
`C0, г/л` = C_0,
`Vd, л/кг` = V_d,
`AUC0-t, г*мин/л` = AUC,
`Cmax, г/л` = Cmax,
`Tmax, мин` = Tmax)
pharmacokinetics_data %>%
pivot_longer(!subject_id,
names_to = "Показатель",
values_to = "Значение") %>%
group_by(`Показатель`) %>%
summarise(across(`Значение`, stat_ph, .names = "{.fn}")) %>%
ungroup() %>%
flextable() %>%
theme_box %>%
fontsize(size=9) %>%
add_header_row(values = "Описательные статистики фармакокинетических показателей", colwidths = 5)
Описательные статистики фармакокинетических показателей | ||||
|---|---|---|---|---|
Показатель | Среднее | Медиана | Стандартное отклонение | Коэффициент вариации (CV), % |
AUC0-t, г*мин/л | 212.3156 | 210.4500 | 25.1503 | 11.85 |
C0, г/л | 0.9886 | 0.9806 | 0.0747 | 7.56 |
Cmax, г/л | 0.9250 | 0.9300 | 0.1550 | 16.76 |
Tmax, мин | 54.3750 | 60.0000 | 28.1263 | 51.73 |
Vd, л/кг | 0.6916 | 0.6935 | 0.0512 | 7.41 |
β, г/л в час | 0.1238 | 0.1243 | 0.0123 | 9.94 |
hist_fn <- function(data, col) {
ggplot(data, aes(x = {{col}})) +
geom_boxplot(
alpha = 0.5,
# width = width,
linewidth = 1.5,
outlier.size = 2
) +
geom_histogram(
fill = "midnightblue",
alpha = 0.5,
colour = "black",
bins = 25
) +
geom_vline(aes(xintercept = mean({{col}}, na.rm = T)),
color = "darkred",
linetype = 2,
linewidth = 2) +
labs(# title = rlang::englue("{{col}}"),
y = "Количество, n",
x = rlang::englue("{{col}}")) +
theme(plot.title = element_text(hjust = 0.5),
axis.title = element_text(size = 25),
axis.text.y = element_text(size = 23),
axis.text.x = element_text(size = 23),
legend.title = element_text(size = 23),
legend.text = element_text(size = 23))
}
col_list <- pharmacokinetics_data %>%
select(-subject_id) %>%
names() %>%
syms()
hist_plots <- purrr::map(col_list,
~ hist_fn(pharmacokinetics_data, !!.x))
wrap_plots(hist_plots, ncol = 2) +
plot_annotation(title = "Гистограммы и боксплоты фарм. показателей",
theme = standard_theme)
Фармакокинетические кривые в целом лежат достаточно кучно, имеют очень похожие профили;
В основном \(Cmax\) достигалась на 30 минуте, всего представлены 4 различных \(Tmax\);
Похоже, что распределения константы элиминации и площади под кривой бимодальны;
Средняя \(Cmax\) составила 0.93 ⼠ 0.16 г/л;
Не смотря на гомогенность выборки по основным параметрам у \(Tmax\) значительный коэффициент вариации (около 50%), возможно, не учтён какой-либо фактор;
Мне кажется, что полученное распределения \(C_0\) и \(β\) имеют достаточно низкую вариабельность.